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Abstract 

A study of the transport coefficients of a system of elastic hard disks, based on the use of Helfand- 
Einstein expressions is reported. The self-diffusion, the viscosity, and the heat conductivity are 
examined with averaging techniques especially appropriate for the use in event-driven molecular 
dynamics algorithms with periodic boundary conditions. The density and size dependence of the 
results is analyzed, and comparison with the predictions from Enskog's theory is carried out. In 
particular, the behavior of the transport coefficients in the vicinity of the fluid-solid transition is 
investigated and a striking power law divergence of the viscosity in this region is obtained, while all 
other examined transport coefficients show a drop in that density range- in relation to the Enskog 
prediction. 

PACS numbers: 05.60.-k,02.70.Ns,05.20.Dd 
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I. INTRODUCTION 



Transport coefficients characterize the different dissipation mechanisms in non- 
equilibrium states. At the macroscopic level, they are introduced by phenomenological 
equations, like the Navier-Stokes equations for a simple fluid, which predict the time evolu- 
tion of mass, momentum and energy [J. Each transport coefficient is related to the propa- 
gation of one (or more) of these microscopic quantities, bridging therefore the hydrodynamic 
and the microscopic scale. In the case of low density gases, the macroscopic equations have 
been justified, their range of validity has been determined, and explicit expressions for the 
transport coefficients have been obtained using the Boltzmann kinetic equation as starting 
point . At higher but moderate densities, the Enskog equation has also proved to give 
a quite accurate description of a gas of hard spheres or disks. 

A remarkable and fundamental development in the statistical mechanics theory of trans- 
port processes was the derivation of expressions for the transport coefficients based on equi- 
librium time- correlation functions. These are the so-called Green-Kubo formulas, and they 

n 

involve different microscopic fluxes s4]. These expressions, although formal, are of general 
validity and have been extensively used for the analysis and modelling of transport in dense 
systems. In particular, they have been applied to the computation of transport coefficients 
by means of molecular dynamics simulations. 

Alternative formal expressions for the transport coefficients are provided by the Einstein- 
Helfand formulas |f|. These are analogous to Einstein's formula for the self-diffusion 
coefficient in terms of the second moment of the displacements. The Einstein-Helfand ex- 
pressions for the other transport coefficients involve moments of dynamical variables which 
are the time integrals of the microscopic fluxes appearing in the Green-Kubo relations. 

In the last years, there has been a revived interest on transport processes in systems com- 
posed by hard particles motivated by the study of granular media in general, and granular 
Erases as a special case [3, 1$, li| The simplest model for them is an assembly of inelastic 
hard spheres or disks, in which the inelasticity is accounted for only through one constant 
parameter, the coefficient of normal restitution. In the low density limit, hydrodynamic 
Navier-Stokes like equations, with explicit expressions for the transport coefficients, have 
been obtained for this model, by starting from the inelastic extension of the Boltzmann 
equation Q|. Moreover, it has been shown that the transport coefficients for a dilute gran- 
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ular gas can be expressed in the form of generalized Green-Kubo relations 
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141 ]. but its density 



The Enskog equation has also been extended to inelastic particles [13J, 
and inelasticity range of validity is not clear. On the other hand, formally exact relations 
between transport coefficients and appropriate correlation functions, similar to the Green- 



Kubo and Einstein-Helfand formulas, appear to be limited up to now to the low density 



limit mentioned above and to the simplest cases of tagged particle motion Jl5l. Il6j]. although 



there have been some more general proposals 0, Il8|| . Therefore, for high densities, the 
only available hydrodynamic theory for granular systems is restricted to the so-called quasi- 
elastic limit. The transport coefficients are given in that limit by the same expressions as 
for elastic systems, and dissipation in collisions is taken into account only by a new term in 
the energy balance equation [toI ]. 

The first calculations of transport coefficients for hard-sphere systems by means of equi- 
ibrium molecular dynamics simulations go back to the pioneer works by Alder and coworkers 
2(1 121I ] . The dependence of the values of the transport coefficients on the density and also 
on the number of particles used in the simulations have been analyzed. At high densities, 
significant deviations from the Enskog theo ry pred ictions are observed, especially for the 
self-diffusion and shear viscosity coefficients |2i 



23, 



24, 



In spite of all the work done for three dimensional systems, it is hard to find results 
for two dimensions, i.e., for a fluid of hard disks. It could be argued that this is due to 
the presence of long time tails in the equilibrium correlation functions appearing in the 
Green-Kubo expressions of the transport coefficients, but it must be noticed that they do 
not invalidate by themselves the possibility of a hydrodynamic description j^fj]. A notable 
exception is ref. (27], where the viscosity of a system of hard disks is measured by using an 
Einstein-Helfand expression. 

Since, in granular materials, a wide range of densities is typically realized, it is useful to 
measure the transport coefficients of a system of hard disks in the whole range of densities. 
In this context, a global equation of state for a system of hard disks has been proposed 28]. 
The equation has proven to be accurate for densities ranging from the low density region 
to the highest crystallization limit where caging effects appear and free volume theories 
are relevant. More precisely, it appears to be almost exact for most densities, except for 
those when the system changes from a disordered state to an ordered one. In this transition 
regime, memory effects become important, hysteresis shows up, and the proposed equation 
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of state is only an approximate description for very slow changes in the density. 

In this paper, the transport coefficients of a system of hard disks will be measured by 
means of Einstein-Helfand expressions that are appropriate for molecular dynamics simula- 
tions with periodic boundary conditions, using the minimum image convention [6] . For con- 
tinuous interaction potentials, this method is strictly equivalent to the Green-Kubo method 
and has the advantages of directly showing the positivity of the transport coefficients and 
of being based on a straightforward, numerically robust accumulation [27]. In the case of 
hard particles, there is a fundamental reason to employ methods based on Helfand-Einstein 
expressions for the transport coefficients. The Green-Kubo relations, except in the case of 
self- diffusion, involve forces between the particles, which are ill-defined for hard spheres or 
disks and there is no trivial way to extend them to hard-particle fluids. In fact, a recent 
careful analysis of the dynamics of a system of hard spheres has shown that the correct 
Green-Kubo expressions for this system have a new singular contribution due to instanta- 
neous collisions, as well as the usual time integral of the flux correlation functions [3] . The 
singular part vanishes in the low density limit, but gives a relevant contribution at high 
densities. On the other hand, the Einstein-Helfand formulas do not involve the forces, and 
have the same form for both continuous (soft) and rigid (hard-sphere) potentials. It must 
be stressed that it is an Einstein-Helfand method that was already numerically implemented 
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by Alder et al. in their study of the transport coefficients in hard-sphere fluids [2( 

The most interesting finding of the present study is that the shear viscosity shows a diver- 
gence at the crystallization density (in a non-sheared system), while the heat conductivity is 
correlated to the isotropic pressure. Thus, while pressure and heat conductivity show a small 
drop at crystallization (due to the better ordering of the particles) , self-diffusion vanishes at 
it, and shear viscosity diverges. This divergence is well below the excluded volume caused 
divergence of pressure and heat conductivity at the maximum density possible in hard disk 
systems. 

The outline of the paper is as follows. In the next Section, the Einstein-Helfand expres- 
sions for the transport coefficients are revised, and written in a way that is appropriate 
for hard sphere molecular dynamics simulations with periodic boundary conditions. The 
special case of self-diffusion, in which the actual trajectories of the particles along different 
cells must be followed, is first discussed. For all the other transport coefficients, it is shown 
that a decomposition of the contributions to the transport coefficients into a kinetic and a 
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collisional part, allows the use of the minimum image convention. Moreover, the decompo- 
sition is especially useful for event driven simulation algorithms. In Sec. IIII1 the method is 
applied to a system of hard disks for calculating the self-diffusion, shear viscosity, and heat 
conductivity coefficients over the whole density range. The results are compared with the 
theoretical predictions from Enskog's theory. Particular emphasis is put on the behavior 
of the transport coefficients in the fluid-solid transition region and on characterizing the 
divergence of the shear viscosity. The paper finishes with a brief discussion of the results in 
Sec. EH 



II. EINSTEIN-HELFAND EXPRESSIONS FOR THE TRANSPORT COEFFI- 
CIENTS 



A. The self-diffusion coefficient 



Self-diffusion is the macroscopic transport phenomenon describing the motion of tagged 
particles in a fluid at equilibrium, in the limiting case that their concentration is very low, 
while at the same time they are mechanically identical to the fluid particles. The macroscopic 
number density n(r, t) and the flux j(r, t) of tagged particles satisfy the continuity equation 

d t n(r,t)+V-j(r,t) = 0. (1) 
The corresponding constitutive relation closing the above equation is provided by Fick's law 

j(r,t) = -DVn(r,t), (2) 

which defines the self-diffusion coefficient D. An expression for this transport coefficient 
in terms of the second moment of the displacements is given by the well-known Einstein 
formula 

where r(t) is the position of an arbitrary tagged particle at time t, d is the dimensionality 
of the system, and the angular brackets denote an average over the ensemble describing the 
equilibrium of the system. 

To actually compute the above equation in our numerical algorithm, two different averages 
have been carried out. First, an average over the N particles in the system is taken, and 



then a second average over a number M of initial configurations (trajectories). Assuming 
ergodicity of the system, different trajectories can be generated from the same simulation 
run by considering different initial times Therefore, the full average is 

1 M N 

(\r(t) - r(0)\ 2 ) m J2 E M + **) " r *(**)| 2 - ( 4 ) 

k=l i=l 

This double averaging is possible because the dynamical variable involved in Eq. © is a 
mono-particle property in the present case, and all the particles are equivalent. 

When using periodic boundary conditions to evaluate Eq. (£Q), as in the simulations to be 
reported here, it is crucial to take into account that particles originally in the center cell in 
a given trajectory, have to be followed as they cross the border of the cell. The positions in 
Eq. (HJ), as writen down here, use the actual positions, but the real relative displacements 
should be used instead. If periodic images of the center cell were not used, the displacements 
would not be obtained correctly. In this sense, the practical implementation of the algorithm 
for this transport coefficient differs from those used for the other transport coefficients to be 
discussed in the following. 

The simulation results for D to be reported later on, will be scaled with the value obtained 
from the Enskog equation in the first Sonine approximation for d = 2 that is j^ : 

De = !-pr (— ) ■ (5) 

2nag2{cr) \ nm J 

Here T is the temperature, ks the Boltzmann constant, a the diameter of the disks, m their 
mass, and (72(c) the value of the equilibrium pair correlation function at contact, which 
is a function of the density n. An estimate for this quantity is provided by Henderson's 
expression |3l|: 

with v = mT<y 2 /A being the solid fraction. The approximation in (JBJ) is valid for densities 
well below the crystallization solid fraction v c ~ 0.7. 



B. Shear viscosity 

The coefficients of shear viscosity 77 and bulk viscosity ( are defined through the macro- 
scopic expression for the pressure tensor P for a simple fluid []] 

P(r, t)=pl-ri [Vu + (Vti) + ] + (^j - C j 1 V • u, (7) 



where p is the pressure, u the flow velocity, 1 the unit tensor, and the superscript + means 
here transposed. 

The Einstein-Helfand formulas for r] and ( are analogous to Eq. (J3J). For the shear 
viscosity one has ^, (| 

with V the volume of the system, and 

N 

G 11 = ^2x i y i . (9) 

i=i 

In this case, only an average over a set of different initial conditions can be taken from 
the simulations, since the dynamical variable G v already involves all the particles in the 
system. Therefore, the quantity that, in principle, should be obtained from the simulations 
is 

N ( N y 

(\G n (t) - a,(o)i 2 > = jr?J2\J2 + tk ^y^ - \ ■ w 

k=i i i=i j 

Nevertheless, this expression leads to very noisy results and the slope of the best fit line 
can be only determined with a high uncertainty. Moreover, it presents the same difficulties 
as the Einstein expression for the self-diffusion coefficient when using periodic boundary 
conditions, i.e., for the positions it requires to follow the motion of the particles through 
different unit cells. It is thus convenient to elaborate a little more Eq. (jlOj) . The idea is 
measuring the increments of G v for physically well defined time intervals, instead of directly 
determining its actual value at successive times. 

Let us consider a time interval [t + tk, t + 1\~ + At] in which no collisions occur in the 
system in a given trajectory k. The purely kinetic variation of G v , due to the displacements 
only, is 

JV 

AGf\t + t k ) = ^2xi(t + t k )yi{t + t k )At. (11) 

i=l 

In addition, there is also a contribution due to the discontinuous change of the velocities in 
collisions. Consider a collision between particles % and j. There is an instantaneous jump in 
G v given by 

AGjpxfyi + x+yj - xTy { - xjy j: (12) 
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where we have taken into account that the positions do not change during the instan- 
taneous collision and the index + (— ) indicates that the velocity is the post-collisional 
(pre-collisional) one. Both velocities are related by the collision rule for hard disks 



= r j + fij ■ era. (13) 

with Tij = Ti — rj and a being the unit vector joining the centers of disks i and j at contact 
and pointing from particle j to particle i. Using the above rule, Eq. (|T2j) can be rewritten 

as 

AGf > = y^, (14) 

where yij = y% — yj, and 8fi = —(rij ■ &)& is the change of the velocity of particle i in 
the collision. Since the dynamics of a hard particle system consists of free streaming and 
instantaneous collisions, Eqs. (jllj) and (|14j) fully determine the time evolution of G n along 
a trajectory. Moreover, these equations only involve the velocities and relative positions of 
the particles. As a consequence, they avoid the difficulties of using other Helfand-Einstein 
relations with periodic boundary conditions in the simulations, as discussed by Erpenbeck 
since no contribution leads to the growth in time of the dynamical variable G v due 
to the infinite checkerboard of identical systems. This is because contributions from pairs 
of particles in different unit cells cancel out precisely due to the boundary conditions. An 
alternative use of a Helfand-Einstein relation for computing the shear viscosity with periodic 
boundary conditions has been discussed in 27 1. 

Equations (fTTj) and ()14|) are particularly suitable for event driven algorithms as the one 
employed in the simulations presented in this paper. In these algorithms, the time steps are 
the intervals between successive hard collisions in the system. At every collision, the kinetic 
change AGq associated with the previous time step is computed as well as the contribution 

(C) 

from the collision itself, AG^ . For the latter, only the positions and velocities of the pair of 
colliding particles must be taken into account, while for the kinetic contribution the motion 
of all the particles in the system has to be considered. 

As for the self-diffusion coefficient, the simulation results for rj will be reported scaled 
with the Enskog value in the first Sonine approximation, that for d = 2 and densities v < v c 
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is 

r i / x\ 1 

(15) 



2v + I 1 + - ) g 2 (a)u 2 



where rj is the value in the Boltzmann limit 

1 / m/c^T\ 1/ 



?7o 



2a V 7T 



J ■ (16) 



C. Bulk viscosity 



The coefficient of bulk viscosityC was already introduced in Eq. (J7J). The corresponding 
Einstein-Helfand expression is P.]^: 

where 

AT 



i=l 



The pVt term in Eq. (|T7j) arises from the fact that the equilibrium average of does not 
vanish, as it is the case for all the other variables G associated to transport coefficients, but 
it is e q na, to the external - P V, deEned by the vhial theo.en, fl. Sinee this m ean 
value is also computed in the simulations and it slightly shifts as the simulations proceed, 
the results for the right hand side of Eq. (|17|) are determined much less accurately than for 
the expressions for the other transport coefficients. Additionally, the subtraction of rj, which 
itself is determined with some uncertainty, causes further errors in the values estimated for 
r). For these reasons, we have not been able to obtain reliable results for the bulk viscosity in 
the high density region, and no further consideration will be given to it here. Note however, 
that it is still possible to split expression (j!7|) in a kinetic and a collisional contribution, like 
we have done above. For that reason, the microscopic expression of the hydrostatic pressure 
is used 0|: 

mN . . . . ma A . 
V = —^\XiXi) + —^(XijdXi), (19) 

where X „ is the unitary vector pointing from center of particle j to center of particle i, and 
T is the collision frequency f2j. In terms of this equation and following the same reasoning 
as in the previous section, it is straightforward to identify that the increments on G^ in the 
particular case of the bulk viscosity, are related to deviations of the variable with respect to 
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the mean value. This happens for the kinetic part, 

N 

AGf\t + At) = J2 [Xi(t)xi(t) - (xi(t)xi(t))} At, (20) 

8=1 

as well as for the collisional part: 

AGf\t + At) = [xijSii - (xijSxi)] , (21) 

with 5xi being the change of velocity of the particle % in its collision with particle j, as 
already defined above. Note that expressions (j2*Uj) and (J5TJ) are again perfectly compatible 
with the minimum image convection. 

We finally reproduce, for the sake of completeness, the Enskog value of the bulk viscosity 
for hard disks |30| ]: 

s^w wy fi } 

7TCT" \ 7T / 

D. Thermal conductivity 

The coefficient of thermal conductivity A is defined by the Fourier law for the heat flux 
q(r,t) [J 

q(r,t) = -XVT(r,t), (23) 
and the Einstein-Helfand expression for it is 

with 

N 



G A = $>^. (25) 



Here is the energy of particle i, 



i=l 



- 2 - (26) 
As it was the case with Eq. (fTUj) . also Eq. (j2*4^) is not appropriate for numerical simulations 

with periodic boundary conditions, since it involves the positions of the particles. Therefore, 

we are going to transform it in a similar way, i.e., measuring separately the increments of the 

dynamical variable G\ between collisions and the collisional contributions. The variation of 

G\ along a trajectory k in a time interval At such that no collisions take place is 



1 N 
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where the isotropy of the equilibrium state has been taken into account. In this way, the 
statistical quality of the simulation results is increased. 

In a collision between particles i and j, the instantaneous change in G\ is 

AG ^ — X^C^ ~\~ XjCj X^e^ ^J^j 

= Xij(ef - e^) = Xij5ei. (28) 

The second equality follows from the total energy conservation in the collisions here consid- 
ered. Using again the isotropy of the system, we can write 

AG (C) = (gg +m ] )5e l (2Q) 

Equations (|2*7)l and ()29|) do not contain absolute positions of the particles, but only their 
relative positions and their velocities. Therefore, we have obtained again an expression that 
does not require to follow the track of the particles all along the simulation, namely to use 
the so-called unfolded particle coordinates. The values of the position and velocities of the 
particles expressed in folded coordinates can be used, without taking care of the cell crossing. 
Enskog's theory prediction for the heat conductivity is 



Ae — Ao 

where Aq is the Boltzmann value 



1 {9 4 , , 

+ 3z/+ - + - irg 2 {<r) 



9i\p) V 4 v 



(30) 



Ao = ^(^Y /2 . (31) 



III. RESULTS 

In this Section, results from several series of event driven simulations are presented. 
We consider an homogeneous, freely evolving system of elastic granular disks with peri- 
odic boundary conditions, implemented by means of the minimum image convention j^]. 
Systems with different numbers of particles and volume have been simulated, and the de- 
pendence of the transport properties on density and particle number has been investigated. 
A typical simulation started with a square lattice of particles having a Gaussian velocity 
distribution. After a transient period, the system reaches a equilibrium homogeneous state. 
Then, the measurement of the different properties of interest was carried out, using the pro- 
cedure described in the previous section. For every system, the above process was repeated 
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a number of times, typically 300, in order to generate the ensemble average over different 
trajectories. Averages over different initial times were also considered for overriding the 
lack of statistical precision in the cases of the shear viscosity and the thermal conductivity, 
as compared with the self-diffusion coefficient. That implied performing longer simulations 
and store and handle the numerical data. The number of initial times used was in all cases 
larger than 200. 



A. Self-diffusion coefficient 

In Fig. [TJ the results obtained for the self-diffusion coefficient, D sim , as a function of 
the solid fraction v are presented. For each density, systems with different numbers of 
particles N have been considered, namely N = 169, 625, 1024, and 2401. Moreover, the 
reported values are the average over 300 independent trajectories. A strong deviation of 
the Enskog value De is observed, even for relatively low densities. Also shown are some 
previous results obtained by Holian et al. jzj sXv = 0.3 using a non-equilibrium molecular 
dynamics method. They consistently agree with the results being reported here. For large 
enough packing fractions z/, the particles become trapped in a crystalline lattice and no free 
movement is possible j^j]. Therefore, the self-diffusion coefficient must vanish in this limit. 
This explains the rather fast decay to zero observed in the simulations. Of course, these 
high density effects are not captured by the Enskog equation. 

The series of values of D obtained for each value of N in the interval < v < 0.5 have 
been fitted to a third degree polynomial 



a + bv + cv l + dir. (32) 



D E {v) 

The values of the fitting coefficients are given in Tabled There, it is seen that the coefficient 
a, characterizing the dilute limit, seems to be weakly dependent on the number of particles 
used in the simulation, and clearly larger than unity, indicating that the dilute limit is slightly 
underestimated by the expression for Dg we have used. Moreover, we have carried out 
simulations with different boundary conditions and found always the same value consistently. 



In fact, similar deviations have been previously observed 36]. This can be due to the 
fact that the Enskog expression given by Eq. (j3J) has been computed in the first Sonine 
approximation, as already mentioned. It is possible that the consideration of higher order 
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FIG. 1: Scaled self-diffusion coefficient for a hard-disk fluid, as a function of the packing fraction 
v. For each value of the density, systems of N = 2401, N = 1024, N = 625, and N = 169 
identical particles have been considered. Also included are previous results obtained by Holian et 
al. using non-equilibrium simulation techniques. Note that the D s i m = De for v = 0.69. The 
extrapolated tangent of the curve at this point crosses the line D s i m = at v ~ 0.705, a value 
which is fairly close to v c . 

polynomial corrections in the Sonine expansion would improve the agreement between theory 
and simulations in the low density limit. 

Figure and Table H] clearly indicate a strong dependence of the measured value of the 
self-diffusion coefficient on the system size. Of course, it is expected that the simulation 
results converge to a well defined value as the number of particles increases, although this 
is not at all clear from Fig. ^ especially for densities around v = 0.5. To check this 
convergence and characterize it, two series of simulations corresponding to v — 0.001 and 
v = 0.5, respectively, have been performed. The results, as functions of the number of 
particles used, are presented in Fig. El It is seen that in the low density case, accurate 
size-independent results are already obtained with a small number of particles, namely with 
iV = 169. On the other hand, for v = 0.5 the convergence is much slower, and reliable results 
require a few thousands of particles. More precisely, the dependence on N in this cases is 
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N 


a 


b 


c 


d 


169 


1.0207 


0.0433 


8.629 


-11.43 


625 


1.0299 


0.3367 


9.775 


-12.39 


1024 


1.0259 


0.4683 


10.711 


-13.48 


2401 


1.0287 


0.5971 


11.930 


-14.75 



TABLE I: Empirical fit of the simulation results for the self-diffusion coefficient to the third order 
polynomial in Eq. (|32|) . The error of the fitting for each value is of the order of the last figure 
given. 
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FIG. 2: Dependence of D S i m on the solid fraction v for low and moderate densities. The symbols 
are some of the simulation results given in Fig. ^ The lines are the fits to the polynomial (|32|) 
with the coefficients given in Tabled In the inset, the low density region is enlarged. 

quite well fitted by an exponential function — D exp(— N/Nq), with N = 900. In any 
case, the existence of an asymptotic value of follows clearly from the above results. 
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FIG. 3: Dependence of the measured self-diffusion coefficient on the number of particles N for 
systems with densities v = 0.001 and v = 0.5, respectively. The results have been averaged over 
30 trajectories for the system with N = 10000 particles, and over 300 trajectories for all the other 
systems. The error bars in this scale are smaller than the symbols used. 

B. Shear viscosity 

The simulation results for the shear viscosity at low and intermediate densities are shown 
in Fig. 0] The deviations from Enskog theoretical predictions is typically under 10% of the 
value for low densities, and this extends up to the transition to the ordered state. This is true 
even for the simulations with the lowest number of particles. In fact, no strong dependence 
of the measured value of the transport coefficient on N can be inferred from the simulation 
data in this range. A similar behavior was found by Alder et al. for a system of hard spheres 

Q. 

Enskog theory clearly underpredicts shear viscosity in the range 0.5 < v < 0.68, contrary 
to what was found in the case of the self- diffusion coefficient, which dropped at v c . This is so 
because the collision frequency in that range of densities is overestimated by Enskog's theory. 
In the crystalline region, the measurement of the shear viscosity becomes rather difficult, 
since the expected linear behavior of the increment in time of the dynamical variable G v (t) 
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FIG. 4: Shear viscosity coefficient normalized with the Enskog prediction as a function of solid 
fraction for dilute and moderately dense systems. The symbols are from the simulations with 
different number of particles N. Averages over 300 different initial conditions and over 300 different 
initial times have been taken. The low density region is amplified in the insert. 

disappears (see Eq. (JHJ). 

In Fig. 03 the deviations of the measured values of T] s i m from Enskog theory in the high 
density region are plotted in a logarithmic scale. A power-law divergent behavior of the 
viscosity is observed as the density approaches the critical value. Moreover, no shift of 
critical (viscosity) density u v is observed as the number of particles increases. The dashed 
line in the figure is the function 



(u) = c(u v - u) 1 , 



(33) 



with c = 0.037±0.001, and u v = 0.71±0.01. This latter value approximately agrees with the 
density for which the self-diffusion coefficient vanishes (see Fig. [TJ, and also with the critical 
(crystallization) density in the global equation of state proposed in Let us mention 

that for = 169 no linear behavior of G n was found for densities v > 0.65 and, therefore, 
no results with this number of particles are included. The behavior of r] sim in this region, in 
fact, seems to depend on N, as can be observed in the figure. 
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FIG. 5: Divergent behavior of the shear viscosity coefficient. The symbols show simulations with 
different number of particles, as indicated. They have been obtained by averaging as in Fig^I The 
dashed line is the power-law fitting given in Eq. (|«S3j) . 

It is worth stressing here the relevance of making the measurements of the corresponding 
dynamical variable in the simulations with a frequency higher than the collision frequency 
in the fluid. In this case, the validity of the numerical procedures discussed in the previous 
section is guaranteed and the results obtained from them can be expected to be correct. 
On the other hand, if the time interval between successive measurements is increased with 
respect to the mean collision time of the system, the results for the time evolution of the 
corresponding dynamical variable G, may not present a linear behavior. But, even if it turns 
out to be linear, the slope can lead to wrong values of the corresponding transport coefficient. 
As an example, a comparison of the measurement of G v in two simulations made with the 
same system of N = 1024 particles and density v = 0.30 is presented in Figure El More 
accurate results are obtained with a time step between measurements of G v shorter than 
the mean time between collisions. This is a consequence of the effect of multiple collisions 
occurring between successive measurements of G v , which invalidates the arguments leading 
to Eqs. (jlljl and (|14j) . Although this applies, in principle, to both the shear viscosity and 
the thermal conductivity, the simulation results show that the influence of the time step 
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FIG. 6: Time evolution of the dynamical variable (AG^) 2 obtained with two different time intervals 
At, one shorter than the mean collision time tE and the other larger than it. The number of particles 
is N = 1024 and the density v = 0.30, in both cases. The solid line is the Enskog prediction. Time 
is scaled by to, Boltzman's mean time between collisions. Note that to = tE in the dilute limit, 
where 52(0") — ► 1- 

between the measurements employed is stronger for the former than for the latter. 
C. Thermal conductivity 

Fig. [7| depicts the results obtained for the heat conductivity. Although the discrepancies 
with the Enskog predictions are not large (they are of the order of 10%), the qualitative 
behavior as a function of the density resembles that of the self- diffusion coefficient reported 
in Fig. ^ It exhibits a maximum around v ~ 0.55, decaying below the Enskog prediction 
for larger densities {y > 0.7). This decay is due to the decreased mean free path due to 
the ordering of the particles. Of course, in contrast with D sim , X sim does not vanish in 
the ordered region, since there is still considerable transport of energy through collisions. 
Moreover, G\ was found to exhibit linear behavior in the transition to the ordered state, even 
for the smallest system considered (N = 169). Similarly to the case of the shear viscosity, 
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FIG. 7: Thermal conductivity coefficient as a function of the solid fraction v for an elastic hard-disk 
systems of different numbers of particles N, as indicated. The values have been averaged in the 
same way as in Fig. 0] Error-bars are shown when larger than the symbols used. 

no systematic dependence of the results on the number of particles used is observed. Let 
us mention that, although small, the deviations from the Enskog predictions in Fig. 0are 
larger than those found by Alder et al. for a system of hard spheres 20]. 

In Fig. [HI we investigate the transition region [v > 0.65), where the dispersion of our 
measurements is clearly higher. The deviation of the measurements with respect to Enskog's 
value is plotted using two different expressions for the pair correlation function ^((x) m 
equation p0]l. On the one hand, we have used the formula given in equation ©, and these 
are the open symbols. For the other two series (solid symbols), a semi-empirical formula as 
given in reference j^, valid in the whole range of densities studied here, was used: 

Pdenseiy^) 



9 2 (f) = 5-2(0-) + m(v\v c , mo) 



2u 



(34) 



where m{v\v c , m <\) is a connecting function, and Pdense is the reduced pressure in the dense 
region (see Ref. |28( for more details). The explicit expressions of these functions are next 
given in terms of the parameters u c , u v and m w 0.0111: 

1 



m(v\v c ,m ) 



1 + exp [-{v - u c )/m ] ' 
19 
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FIG. 8: Deviation of the values measured for the thermal conductivity coefficient, A s i m , from 
Enskog's value, A^, as a function of the solid fraction v for dense systems. The values have been 
averaged in the same way as in Fig. HJ For the scaling of the tagged series N = 1024* and 
N = 2401*, the values v c = 0.71 and mo = 0.01 have been used. More details about the scaling of 
these results are given in the text. The solid line indicates the fixed value (Xsim/^E — 1) = 0.15. 



Pdense{v) ~ [l - 0.04(i/„ - u) + 3.25(i/„ - uf] . (36) 



v n -v 



In Fig. |H1 results for high values of v are plotted. The log-log representation is used in 
order to be consistent with Figure El The empirical pair correlation function g% is expected 
to work better for values of the density above v w v c = 0.70 [3,0,12^. The deviation from 
the theoretical value remains approximately constant for a wide range of solid fractions, 
covering the low and moderate densities. When g% is used, the deviation remains constant 
even beyond the transition to the dense configuration, while there is a clear deviation of 
the data if the simplified form of equation (JOJ) is used. This is more clearly observed in the 
range (0.65 < v < 0.74). The effect of the correction for smaller values of v is nevertheless 
negligible. 



It has been suggested 



23| that the good quantitative agreement between the simulation 



results and the Enskog prediction for the heat conductivity, may be attributed to the absence 
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of a significant "tail" in the autocorrelation function appearing in its Green-Kubo expression. 
Nevertheless, the Helfand-Einstein approach of the transport coefficients we have used, is 
free from the practical difficulties associated with those tails. Moreover, the fact that a 
linear time dependence of the relevant Helfand moment has been found for both the shear 
viscosity and the heat conductivity, seems to confirm the validity of the formulas used. 



IV. DISCUSSION 



In this paper, the transport coefficients of a system of elastic hard disks have been eval- 
uated by means of equilibrium molecular dynamics simulations. In order to avoid the fun- 
damental difficulties rece ntly identified in the use of the standard Green-Kubo formulas 
in hard-particle systems 29(, a Helfand-Einstein representation has been employed. The 
consideration of periodic boundary conditions in the simulation of a non-sheared, isotropic, 
homogeneous, freely evolving system of elastic hard spheres or disks forces some modifica- 
tions of the usual Einstein- Helfand's formulas for the transport coefficients, especially if the 
minimum image convention is used. Moreover, the expressions proposed here are especially 
suitable for event driven methods. They allow a detailed study of the dependence of the 
coefficients on the system size and density. 

For the self-diffusion coefficient, D, the Enskog approximation leads to values that un- 
derestimate the simulation results by factors up to two, for moderate values of the density 
[y < 0.3), the discrepancies being already relevant at rather low densities. The observed 
density dependence of the transport coefficient is well fitted by a third order polynomial for 
v < 0.5, with coefficients that slightly depend on the number of particles, N, of the system. 
It has been verified that D tends to a well defined limit as iV becomes large enough. At 
higher densities, the transition liquid-solid is clearly depicted in the behavior of the self- 
diffusion coefficient. It rapidly falls to zero as a consequence of the caging of the particles. 
Finite size effects are more relevant for dense systems, in which the self-diffusion coefficient 
approaches its asymptotic value exponentially with N. 

For the shear viscosity the dependence of the results on the size of the system is much 
smaller. Also much weaker deviations from the Enskog prediction are observed at low and 
intermediate densities. Nevertheless, closer to the gas-solid transition, a power law divergent 
behavior has been identified. Interestingly, the density value for which the viscosity would 
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FIG. 9: Schematic plot of the transport coefficients. The solid line gives the pressure P = l + 2vgL, 
with gi from Eq. 1)351) ; the dashed line gives the scaled shear viscosity, i.e. Eq. 1)16)1 . where g^ is 
used instead of g<i , and multiplied by our empirical correction factor 1 + in Eq. 1)34)1 ; and the 
dimensionless heat condictivity from Eq.(31), also with gi used instead of g2- 

become eventually singular [y c ~ 0.71), agrees with the density at which the system begins 



to show an ordered triangular structure 
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4l|. 



Note that the results here presented refer to "non-sheared" systems. We have avoided 
therefore the problem of the system becoming inhomogeneous and developing a shearband 
3|. A sheared system will not show the divergence found for the viscosity because of the 
shearbanding instability. 

The pressure also diverges but at a considerably higher density u max m 0.9069. At the 
crystallization density u c ~ v v ~ 0.7, both pressure and heat conductivity show a drop 
relative to the Enskog prediction due to the better ordering in the crystalline phase (see 
Fig. The use of a more elaborated expression of the pair correlation function, valid in 
a wider range of densities than Henderson's approximation, improves the agreement of our 
data with Enskog theory. There were no data obtained for the shear viscosity above u v . 

In conclusion, we have found Enskog theory working rather well for pressure, heat- 
conductivity and shear viscosity well below the crystallization density v c . When the Enskog 
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expressions are corrected by an appropriate pair correlation function gi, which accounts 
for the better ordering in the crystalline phase, the theory performs well for pressure and 
heat-conductivity up to the maximal possible density v max . 

Only the shear viscosity shows a power-law divergence at v n ~ v c with values above 
Enskog theory already becoming visible at intermediate densities. Thus, shear viscosity 
behaves differently than the other transport coefficients studied. Its divergence, implying 
that the shear modes are hindered for v > u v . This could in fact be understood as one 
reason for shear-band formation. A sheared system at high densities typically splits into 
shear bands (with lower density) and a compressed dense crystal (with correspondingly 
higher density). From a different point of view, our observations are also consistent with the 
concept of dilatancy: A dense packing with v > v v can only be sheared by first experiencing 
dilatancy so that u drops below u c . 
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